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SYSTEMATIC  SPLITTING  OF  WAVEFIELDS  INTO  UNIDIRECTIONAL 
MODES:  LONG-RANGE  ASYMPTOTIC  METHODS  FOR  WEAKLY 

NONUNIFORM  MEDIA 


1.  INTRODUCTION 

This  work  aims  to  produce  a  unified  understanding  of  the  three  analytic  approximations 
WKB,  parabolic,  and  Born  —  that  are  most  widely  cited  to  justify  approximating  wave  propagation 
through  media  with  only  “weak”  nonuniformities  as  a  purely  one-way  phenomenon.  The  particular 
focus  is  on  backscatter  —  how  to  neglect  it  in  a  consistent  way,  and  what  impact  this  neglect 
has  on  the  remaining  forward-propagating  component.  To  limit  the  complexity  of  the  problem 
and  yet  retain  most  of  the  essential  features,  attention  is  restricted  to  waves  in  a  single  spatial 
dimension.  The  subject  is  introduced  in  section  2  in  the  context  of  the  prototypical  example 
vibrations  of  a  nonuniform  string.  This  section  explores  the  difficulties  inherent  in  any  attempt  to 
generalize  forward/backward  mode  splitting  beyond  uniform  media,  where  it  is  a  trivial  matter,  to 
nonuniform  media,  where  it  can  generally  be  done  only  approximately.  In  section  3,  the  initial- 
value  problem  for  the  Helmholtz  equation  is  first  recast  in  first-order  form  and  then  transformed 
by  a  rotation  into  the  “d’Alembert”  representation  where  questions  of  mode  coupling  are  more 
naturally  addressed.  Section  4  discusses  the  effect  of  such  state-space  transforms  on  the  equation 
of  motion  —  especially  the  class  of  “pseudo-unitary”  rotations  that  strictly  conserve  the  wave 
energy.  A  pseudo-unitary  rotation  is  used  in  section  5  to  go  to  the  “Bremmer”  representation, 
where  it  is  evident  what  will  be  required  to  decouple  the  counter-propagating  wave  modes  to 
first  order.  Section  6  introduces  a  scheme  of  successive  pseudo-unitary  Pauli-space  rotations  that 
reveals  the  conditions  for  decoupling  the  modes  to  any  desired  order,  giving  operational  meaning 
to  the  phrase  “weakly  nonuniform  medium”.  In  section  7,  this  weak  nonuniformity  is  invoked,  and 
the  resulting  asymptotic  approximations  are  obtained  for  the  long-range  wave  field  through  order 
m  =  6.  The  phase  in  this  expression  is  seen  to  agree  with  the  WKB  result  through  order  m  —  3, 
with  differences  appearing  at  m  =  4.  Endpoint  amplitude  effects  are  also  discussed,  and  it  is  shown 
that  these  begin  to  differ  from  the  WKB  approximation  at  third  order.  Finally,  section  8  presents 
an  illustrative  example  involving  a  specific  environment.  Beginning  with  section  3,  Pauli  matrices 
are  used  to  facilitate  the  analysis.  Appendix  A  reviews  their  properties  —  particularly  the  fact 
that,  when  used  as  the  infinitesimal  generators  of  state-space  transforms,  the  Pauli  matrices  just 
induce  rotations.  The  relation  of  this  work  to  the  Born  series  approach  is  discussed  in  Appendix  B. 

2.  BACKGROUND 

This  section  sets  the  scene.  The  Helmholtz  equation  for  continuous  wave  (cw)  motion  and  its 
ancillary  relations  for  energy  density  and  power  flux  are  obtained  for  one-dimensional  nonuniform 
media,  and  the  basic  requirement  for  approximate  mode  splitting  is  identified. 

Manuscript  approved  November  14,  1995. 
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Consider  the  transverse  motion  of  a  string  stretched  along  the  x  coordinate  under  tension  r. 
Suppose  the  string  has  a  density  p(x)  and  consider  its  displacement  w(x,t)  in  the  transverse  y 
direction.  Within  the  linear-response  regime,  w  will  obey  the  familiar  linear  wave  equation.  If  the 
x  dependence  of  p  is  only  a  gradual  one  relative  to  the  wavelengths  involved,  the  solution  to  be 
expected  on  physical  grounds  should  consist  of  a  pair  of  waveforms  that  travel  in  opposite  directions 
along  the  string,  changing  only  gradually  with  x  and  only  weakly  coupled  to  one  another. 

Uniform  Medium 

If  r  and  p  are  just  constants,  the  displacement  obeys  a  simple  form  of  the  wave  equation 

pdfw  -  rd2xw  =  0  ,  (1) 

where  the  phase  speed  c  =  \frjp  is  constant.  (See  various  elementary  books,  e.g.  Ref.  1.)  The 
shorthand  notation  dqr  =  dr/dq,  dqr  =  dr/ dq  is  used  for  derivatives  throughout  this  work. 

The  kinetic  energy  is  an  inherently  local  quantity  (because  the  mass  is)  with  a  density 
Kw  =  (p/2)(dtw)2.  The  potential  energy  can  also  be  regarded  as  localized,  with  a  density 
T4)  =  (t/2 )(dxw)2.  The  total  energy  density  for  the  wave  motion  is  just  their  sum, 

Ew  =  ~( dtw )2  +  ^( dxw )2  .  (2a) 

The  power  flux,  also  known  as  energy  current  density,  is 

Pw  =  —T(dxw)(dtw)  .  (2b) 

Since  the  medium  is  passive  and  lossless,  wave  energy  is  conserved, 

dxPw  +  dtEw  =  (dtw)(pd2w  —  rdxw)  =  0  .  (3) 

Consider  the  distribution  of  energy  in  a  wave  field  that  is  a  superposition  of  two  others:  w  = 
a+p.  Since  the  total  wave  energy  is  conserved,  one  might  naively  expect  that  Pa+p  and  Ea+p  would 
simply  reduce  to  Pa  +  Pp  and  Ea  +  Ep,  but  this  cannot  be  true  in  general;  linear  superposition  does 
not  apply  to  energy  quantities  because  they  are  quadratic  in  the  wave  field.  Thus,  the  differences 

Ea,p  d=  Ea+p  -  ( Ea  +  Ep)  =  p(dta)(dtp)  +  r(dxa)(dx(3)  (4a) 

Pa,p  d-  Pa+p  ~  (Pa  +  Pp)  =  - T(dxa)(dtP )  -  r(dxp)(dta)  (4b) 

do  not  vanish  identically.  However,  there  is  one  important  case  where  this  energy  superposition 
does  hold.  When  d’Alembert’s  decomposition  into  left/right-going  waveforms  is  used,  i.e., 

w(x,t)  =  a(£)  +p(rj)  ,  (5) 

where  £  =  x  —  ct  and  rj  =  x  +  ct,  then 

Ea  =  r(d^a)2  >0  PQ  =  +cEa  >  0 

Ep  =  rid^p)2  >0  Pp  =  -cEp  <  0 

Ea+P  =  EQ  +  Ep  >  0  Pa+P  =  PQ  +  Pp 


(6) 
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so  that  both  Ea$  and  PQ)/g  vanish. 

It  is  also  well  known  [1]  that  w  satisfies  Eq.  (1)  if  and  only  if  d’Alembert’s  decomposition 
holds,  i.e.,  that  this  decomposition  provides  the  general  solution  for  the  wave  equation  in  a  uniform 
medium.  Any  wave  field  in  such  a  medium  is  simply  a  pair  of  counter-propagating  waveforms, 
the  d’Alembert  modes,  and  these  partition  the  total  energy.  Energy  superposition  applies,  and 
each  mode  separately  transports  its  share  of  the  total.  These  modes  move  at  constant  speed, 
propagate  without  distortion,  and  partition  the  energy  exactly.  Because  there  is  no  environmental 
inhomogeneity  to  couple  the  modes,  there  can  be  no  backscatter.  If  the  field  is  initialized  with  only 
one  of  these  modes  excited,  the  other  one  can  never  arise  anywhere.  The  aim  of  this  work  is  to 
extend  this  result  in  a  controlled,  approximate  way  to  media  with  weak  inhomogeneities,  providing 
a  generalized  form  of  d’Alembert  decomposition. 

Nonuniform  Media 


The  energy  expressions  in  Eq.  (2)  remain  valid  even  when  r  and  p  are  not  simply  constants. 
The  Lagrangian  density  L(w,  dxw,  dtw )  =  Kw  —  Vw  is  still  governed  by  the  appropriate  Lagrange 
equation 


d 

r  dL  1 

d 

_i _ 

r  dL 

dL 

dt 

d{dtw)_ 

dx 

d(dxw ). 

dw 

(7) 


i.e., 


dt  ( pdtw )  +  dx  (-rdxw)  =  0  .  (8) 

With  p  =  p(x)  and  even  r  =  r(t),  this  would  reduce  to  Eq.  (1)  again.  A  time-dependent  tension, 
however,  would  introduce  a  source  term  on  the  right-hand  side  of  Eq.  (3),  so  we  will  confine  ourselves 
to  cases  with  constant  r.  In  these  static  nonuniform  environments,  wave  energy  is  still  conserved 
—  Eq.  (3)  remains  valid  —  but  d’Alembert’s  decomposition  no  longer  holds  in  the  original  form. 
It  must  be  generalized. 

Complex  Representation 


This  generalization  will  be  easier  to  produce  using  a  complex  representation.  Physically,  the 
wave  field  is  certainly  a  real-valued  quantity.  But  since  the  wave  equation  is  linear  and  has  real 
coefficients,  there  is  no  harm  in  adding  on  some  imaginary  part  to  extend  the  physical  field  wre  to 
complex  values:  w  =  wre  +  Naturally,  both  wre  and  Wim  must  satisfy  Eq.  (1),  and  the  same 

is  true  of  w  and  w*.  Also,  the  dxw,  dtw  factors  in  the  energy  expressions  now  refer  to  dxwre,  dtwre 
so  that 


Ew  - 


|( dtw  -(-  dtw *)2  +  dxw  +  dxw *)2 


/  4, 


Pw  =  -r{dxw  +  dxw*)(dtw  -f  dtw*)/ 4 


(9a) 

(9b) 


As  a  result, 

dxPw  +  dtEw  —  ( dtw  +  dtw*)  ^(pdfW  -  rd^w)  +  (pd^w*  -  rd^w*) J  /4  =  0  ,  (10) 


as  expected. 
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For  the  superposition  field  w  =  a  +  (3, 

Ea,0  =  [p(dta  +  dtam)(dtfi  +  ft/T)  +  r(dxa  +  dxa*)(dx0  +  dx0*)]  /4  ,  (11a) 

PQ, 13  =  [-r(dxa  +  dxa*)(dt(3  +  dtp*)  -  r(dxp  +  dxpT){dta  +  ft  a*)]  /4  ,  (lib) 

but  the  right-hand  sides  fail  to  vanish  because  there  is  no  d’Alembert  decomposition  to  separate  the 
field  exactly  into  counter-propagating  modes.  However,  there  is  a  workable  approximate  approach. 

CW  Fields 


It  might  be  possible  to  deal  successfully  with  fields  whose  time  dependence  is  quite  general. 
Narrowband  fields,  for  example,  could  be  handled  using  analytic-signal  methods.  Since  this  would 
be  too  great  a  diversion  from  the  present  purpose,  however,  the  scope  of  this  work  is  confined  to 
continuous  wave  (cw)  fields. 


Since  the  medium  itself  is  not  time-dependent,  each  frequency  component  can  be  treated 
separately.  With  only  a  single  frequency  u>  in  its  spectrum,  the  wave  field  can  only  be  si¬ 
nusoidal  in  time:  wre(x,t)  =  F(x)  cos(ut  —  <p(x)).  If  the  imaginary  part  is  chosen  to  be 
Wim(x,t)  =  —F(x)  sin(wt  —  <p(x)),  then  the  complex  representation  of  the  wave  field  is 

w(x,t)  =  u(x)e~'ut  ,  (12) 

where  u(x)  —  F(x)e,lfi(x\  Thus  the  wave  equation  for  wre  =  (w  +  w*)/2  reduces  to 

e~iui(u"  +  k2u )  +  e+^V"  +  k2u *)  =  0  ,  (13) 

where  the  primes  denote  x-derivatives  and  k(x)  =  w/c(x),  so  both  u  and  u*  must  satisfy  the 
one-dimensional  Helmholtz  equation 


The  field’s  power  flux  is  then 


where 


and  its  energy  density  is 


where 


u"  +  k2u  =  0  . 

(14) 

Pw  =  Pw  +  Pw  +  Pw  > 

(15) 

P+  =  -e+i2“‘(u*'u>r/ 4 

(16a) 

P~  =  +e~,2u)t  (u1  u)iur  /  4 

(16b) 

P®  =  ( uu *'  —  u*u')iu)T / 4  , 

(16c) 

Ew  =  E+  +  E^  +  Ew  , 

(17) 

E+  =  e+i2ut(u*'2  -  k2u*2)r/ 8 

(18a) 

E~  =  e~i2ut(u'2  -  k2u2)r/ 8 

(18b) 

Ew  =  ( k2u*u  +  u*'u')t/ 4  . 

(18c) 
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Here  we  apply  time  averaging  (denoted  by  a  bar).  This  effectively  kills  off  the  time-dependent  *+’ 
and  ’  terms,  leaving  only  the  constant  ‘0’  terms.  Thus 

Pw  =  W(u,u*)iuT/ 4  (19a) 

Ew  =  ( k2u*u  +  u*'u')t/ 4  ,  (19b) 

where  W  is  a  Wronskian.1  It  is  clear  from  Eq.  (19)  that,  although  Ew  can  depend  on  x,  Pw  cannot. 
That  is  a  statement  of  conservation  of  wave  energy  for  the  cw  case. 


Each  component  of  a  superposition  field  w  =  a  +  (3  can  be  represented  in  complex  form, 
a(x,t)  =  a(x)  exp(— iuit)  and  /3(x,t)  =  b(x)  exp(— iut).  Then, 


pa,0  ~  Pa,/3  +  Pa,p  +  Pa,0  > 

(20) 

where 

P+P  =  -e+i2wt  (a*'b*  +  a*b*')  iur/ 4 

(21a) 

Pa/3  =  +e~,2a,i  (a'b  +  ab ')  icoT/4 

(21b) 

P%  =  [(«*'&  +  -  (a'b*  +  a*6')]  twr/4  , 

(21c) 

and  thus 

Pa,n  =  [W(b,  a*)  +  W(a,  6*)]  iwr/4  . 

(22) 

This  clearly  shows  what  would  be  needed  to  partition  the  power  between  two  modes  a  and  /?. 
a  and  a*  must  be  solutions  to  Eq.  (14),  Pa  cannot  depend  on  x.  The  same  is  true  of  b  and  b*, 
cannot  depend  on  x  either.  Finally,  in  order  for  PQtp  to  vanish  identically,  there  must  be  a 
dependence  between  a  and  b*  (and  thus  between  a*  and  b): 

Since 
SO  Pfj 
linear 

a  oc  6*  . 

(23) 

Similarly, 

E<*,P  -  Etf}  +  E%  +  Ea,p  > 

(24) 

where 

E+p  =  e+i2u>t  (o*'6*'  -  k2a*b*)  r/4 

(25a) 

E~p  =  e~t2uit  {a'b'  —  k2ab^  r/4 

(25b) 

E%  =  [(o'6*'  +  a*'6')  +  k2(a*b  +  ab*)]  r/4  . 

(25c) 

Again,  under  time  averaging,  the  ‘+’  and  ’  terms  vanish,  leaving  only 

EaJJ  =  3?  { a'b *'  +  k2ab*}  r/2  .  (26) 


1Recall  that  W(u,  v)  =  uv  —  vu'  is  independent  of  x  when  u  and  v  are  solutions  to  Eq.  (14),  and  vanishes  whenever 
these  solutions  are  linearly  dependent,  i.e.,  when  roc  b. 
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From  Eq.  (23)  it  is  clear  that  b*  =  where  (  is  some  complex  x-independent  quantity.  Thus, 


At  the  same  time, 


Ea,/3  =  9?  [c(fl/2  +  fc2a2)]  t/2  . 

(27) 

E'a  =  [a'(o*"  +  k2a*)  +  a*'(a"  +  k2a )]  r/4 

(28a) 

E'p  =  [b'(b*"  +  k2b*)  +  b*'{b"  +  k2b )]  r/4  . 

(28b) 

This  is  where  difficulties  appear  in  this  attempt  to  separate  the  field  exactly  into  counter-propagating 
modes  that  simultaneously  obey  the  Helmholtz  equation  and  partition  the  energy.  Because  a,  b, 
a*,  and  b*  are  all  solutions  to  the  Helmholtz  equation,  both  E'a  and  E'p  vanish  identically,  leaving 
Ea  and  Ep  independent  of  x.  This  cannot  be  right;  Ew  itself  did  not  have  to  be  independent  of  x. 
Furthermore,  Ea $  can  vanish  only  when  a'  =  ±ika.  That  implies  a  functional  form  for  a, 


a(x)  =  a(x0)  exp  |±i  J  dxk(x) 
Unfortunately,  it  also  implies  that 


a"  —  k2a  =  ±ik'a  . 


(29) 


(30) 


Since  the  right-hand  side  must  vanish,  it  is  essential  that  k'  =  0.  Clearly,  the  exact  mode  separation 
of  the  field  can  be  done  only  in  a  uniform  medium,  where  it  is  nothing  more  than  the  familiar 
d’Alembert  representation.  It  cannot  be  done  in  nonuniform  media  —  at  least  not  in  an  exact  way. 
Equation  (30),  however,  suggests  that  mode  separation  might  be  achieved  in  some  approximate 
sense  based  on  the  condition 


\k'\/k 2  <  1  .  (31) 

The  remainder  of  this  report  is  directed  toward  realizing  that  possibility. 

3.  INITIAL  REFORMULATION 

This  section  transmutes  the  whole  wave  motion  problem  into  a  form  where  the  effects  of  envi¬ 
ronmental  inhomogeneity  on  mode  splitting  are  clearer  and  easier  to  begin  dealing  with.  Initially, 
the  Helmholtz  equation  and  its  power  flux  and  energy  density  relations  are  recast  in  first-order  form 
involving  Pauli  matrices.  Then  the  Helmholtz  equation’s  Picard  series  solution  is  presented  and 
the  reason  for  its  generally  slow  convergence  is  underscored.  Finally,  the  problem  is  transformed  to 
the  “d’Alembert”  representation  —  an  optimal  starting  point  for  subsequent  perturbation  analysis. 

The  one-dimensional  Helmholtz  equation  has  the  form 

d2u  +  k2u  =  0  , 

where  the  wavenumber  k  —  u/c  may  be  x-dependent.  The  complex  field  is  w(x,  t)  =  u(x)  exp (-iut) 
and  the  relevant  energy  quantities  are 

Pw  =  (— iur/4)  (u*dxu  —  udxu*) 

Ew  =  (t/4)  J k2u*u  +  (dxu*)(dxii) j  . 


Systematic  Splitting  of  Wavefields 


First-Order  Form 

It  is  convenient  to  adopt  some  constant  density  value  po  as  a  standard  and  then  use  the 
corresponding  phase  speed  co  =  y/rjpo  to  define  a  reference  wavenumber  ko  —  u/cq.  The  spatial 
coordinate  will  be  rescaled  to  kox  (so  that  distances  are  measured  out  in  reference  wavelengths), 
but  this  nondimensional  coordinate  will  continue  to  be  written  as  x.  Thus  the  Helmholtz  equation 
becomes 


d2u  4-  v?u  =  0  ,  (32a) 

where  n  —  cq/c  =  k/ko  =  y/p/po  is  the  refractive  index  and,  with  constants  lumped  into  7  = 
a;2po/4,  the  energy  terms  are 

Pw  =  — *co7  ( u*dxu  —  udxu*)  (32b) 

Ew  =  7  [n2tt*u  +  (dxu*)(dxu)]  .  (32c) 


Since  Eq.  (32a)  is  linear,  the  familiar  trick  of  swapping  order  for  dimensionality  [3]  can  be  used 
to  convert  it  from  a  second-order  one-dimensional  form  to  a  first-order  two-dimensional  one.  To 
facilitate  dealing  with  the  complex  2-by-2  matrices  that  the  latter  form  entails,  we  invoke  the  Pauli 
spin  matrices  whose  properties  are  reviewed  in  Appendix  A. 


Initial  Representation 


The  simplest  way  to  produce  a  first-order  form  is  to  let  the  dependent  variables  be  u  and 
u  =  dxu  (the  displacement  of  the  string  and  the  slope  of  its  tangent)  so  that  the  state  of  the  system 
is  specified  by  the  vector 


u  = 


u 

u 


Then  the  Helmholtz  equation  takes  the  form 

dx  u  =  G  u  , 

where 


G  = 


0  1 
— n2  0 


1  4-  n2 


=  i 


o  2  + 


1  —  nr 


<Tl 


(33) 


(34a) 


is  the  generator  of  the  system’s  evolution  along  the  x  coordinate. 


This  is  equivalent  to  the  Schrodinger  equation,  idx u  =  Hu,  for  a  two-state  quantum  system 
with  a  Hamiltonian  H  =  iG  where  x  plays  the  role  of  time  (in  a  system  of  units  where  %  =  1). 
The  Hamiltonian’s  eigenvalues,  ±n,  correspond  to  time-dependent  energy  levels.  The  power  flux 
and  energy  density  are  “matrix  elements”  (in  the  quantum  mechanical  sense) 

P(  u)  =  co7M^P  u 
E(  u)  =  7M^Eu 
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of  a  pair  of  Hermitian  operators, 


P  = 


0 

i 


(34b) 


E  = 


n2  0 

0  1 


(34c) 


The  only  peculiar  aspect  of  this  as  a  quantum  mechanics  problem  is  that,  except  for  the  degenerate 
n  =  1  case,  H  is  non-Hermitian.  The  work  below  is  generally  done  in  terms  of  G,  with  H  appearing 
only  when  analogies  to  quantum  mechanics  are  pursued. 


In  a  certain  sense,  the  solution  is  already  in  hand.  When  Eq.  (33)  is  combined  with  any 
initial  condition  u(xo)  (forming  an  “initial- value”  or  “one-point  boundary  value”  problem),  the 
state  vector  at  any  x  is  simply 


u(z)  =  K(:e,2:o)u(xo) 

in  terms  of  the  propagator  matrix  K(x,x0)-  This  matrix  operator  is  defined  by  a  Dyson-ordered 
integral  [2]  (also  known  as  a  “product  integral”  [4])  and  is  variously  called  a  “matricant”  [5]  or 
“matrizant”  [6].  It  is  the  solution  to  the  same  ordinary  differential  equation  with  unit-matrix 
initial  values: 


dxK(x,  s0)  =  G(s)  K(x,  s0)  ...  K(s0,  s0)  =  1  ,  (35) 

or,  equivalently,  the  integral  equation 

K(x,  so)  =  1  +  f  d£  G(0  K(£,  so)  .  (36) 

Jx  0 

Under  fairly  general  conditions,2  this  has  a  unique  solution  in  the  form  of  an  infinite  series 


K(x,  s0)  =  Y,  K»(*>  s0)  (37a) 

j= 0 


that  begins  with 

Ko(s,s0)  =  1  (37b) 

and  continues  by  Picard  iteration 

Kj+i(x,x0)=  f  d£G(f)Kj(C,*o)  ,  (37c) 

Jxq 

converging  uniformly  and  absolutely  [5,6].  This  Picard  series3  is  satisfying  but  less  useful  than  it 
may  appear,  especially  in  nonuniform  media,  because  the  convergence  is  typically  quite  slow. 

Sufficient  conditions  are  that,  in  an  interval  containing  x  and  xo,  n2(x)  be  single- valued,  bounded,  and  integrable 
and  that  dxn2(x)  be  piecewise  continuous  and  bounded  [6], 

3It  is  also  associated  with  the  names  Liouville  and  Neumann. 
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To  see  why  that  is,  let 

e  =  (n2  -  l)/2  (38) 

measure  the  deviation  of  the  refractive  index  from  its  reference  value,  1,  so  that  the  three  operators 
of  Eq.  (34)  are 


G  =  — e<Ti  +  i(l  +  e)<72 

(39a) 

P  =  <t2 

(39b) 

E  =  (i  +  e)l  +  e<T3  . 

(39c) 

The  generator  G  has  a  constant  term  as  well  as  e-dependent  ones.  Thus  K j,  a  multiple  integral  of 
a  product  of  the  form 


£(*,•) £(*i-0  •••G(*2)G(xi)  , 


contains  contributions  of  multiple  orders  in  e(x).  The  integrand  for  K3,  for  example,  is 


G(x3)G(x2)G(xi) 


0  -(l  +  2e(x2)) 

(1  +  2e(a:3))(l  +  2e(xi))  0 


(40) 


which  has  contributions  of  orders  0,  1,  and  2  in  e.  The  convergence  of  Eq.  (37)  is  slow  because 
Kj  generally  contains  terms  in  e°,  e1,  e2,  •  ■  ■ ,  eJ_1,  eJ.  Although  it  is  mathematically  allowable  to 
expand  all  the  terms  and  re-order  the  series  so  that  the  jth  term  contains  only  eJ,  no  method  is 
known  for  doing  that  with  any  generality.  There  is,  however,  a  way  to  finesse  a  solution  by  applying 
Pauli-space  rotations. 


d’Alembert  Representation 


Before  pursuing  that  option,  it  will  be  convenient  to  change  the  representation  to  one  in  which 
the  generator  G  is  diagonal  wherever  the  inhomogeneity  e  vanishes.  This  is  done  by  applying  the 
unitary  operator 


1  1 
—i  i 


(41) 


to  transform  state  vectors  and  operators  according  to  u  =  Uu  and  V  =  UVU  *,  respectively.  The 
form  of  the  problem  remains  the  same, 


(42a) 
(42b) 
(42c) 

in  terms  of  the  transformed  state  vector 

u  +  iu 
u  —  iii 

and  transformed  operators  G,P,E  in  this  new  “d’Alembert”  representation. 


dx  u  =  Gu 
P(u)  =  co7u*Pu 
E(  u)  =  7U^Eu 
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Equation  (41)  can  be  written  as  U  =  exp{-iS}  in  terms  of  an  operator  S  =  Sol  +  S  <T 
that  has  a  scalar  part  So  =  —  7r/4  and  a  Pauli- vector  part  S  =  Ss  with  magnitude  S  =  tt/ 3  and 
direction  s  =  (ei  —  ei  +  £3)/ a/3.  Thus,  U  transforms  operators  into  the  d’Alembert  representation 
by  simply  rotating  their  Pauli-vector  parts  about  the  (1,  -1, 1)  direction  through  the  angle  -2 S  = 
— 27r/3.  That  rotation  amounts  to  a  cyclic  permutation  of  (e%,  -<23,  — ei),  i.e.,  to  the  replacement 
(ei,e2,e3)  — ►  (-£2,  -e3,  ei),  which  means  that  the  effect  of  U  can  be  found  by  inspection.  In 
particular,  the  three  operators  of  Eq.  (39)  become 


G  —  — i(  1  +  e)(T  3  +  e<T2 

(43a) 

P  =  -<73 

(43b) 

E  =  (1  +  e)l  +  e<Ti  . 

(43c) 

Their  scalar  and  Pauli-vector  parts  are 

Go  =  0  G  —  -z(l  +  e)e3  +  ee>2 

Po  =  0  P  =  -h 

Eq  =  1  +  e  E  =  ee  1  , 

and,  whh  the  positive  sign  chosen  for  square  roots,  the  magnitudes  and  directions  of  the  vector 

parts  G  =  Gg,  P  =  Pp,  and  E  =  Ee  are 

G  —  in  g  = -[(l  +  e)e3  +  iee2]/n 
P  —  1  p  —  S3 

E  —  e  e  =  e\  . 

On  first  inspection,  the  initial-value  problem  seems  no  better  behaved  than  before,  since  the 
generator  in  Eq.  (43a)  will  also  produce  e-ordering  problems  in  the  Picard  series  for  its  propagator. 
The  advantage  of  this  representation  is  that,  within  any  i-interval  where  k  =  fc0,  the  generator 
reduces  to  a  diagonal  form:  G  — >  —i<T3-  Since  the  forthcoming  solution  will  emerge  from  a 
perturbation  process  for  |e|  <C  1,  it  is  naturally  simpler  to  begin  in  the  d’Alembert  representation 
where  the  unperturbed  e  =  0  problem  is  diagonal. 

Mode  Separation 

At  this  point,  it  is  worth  noting  that  wherever  k  =  &o,  the  d’Alembert  generator’s  eigenvectors 
are  |<?)  with  corresponding  eigenvalues  -gi,  where  g  =  ±.  This  means  that  the  state  vector  separates 
exactly  into  a  sum  of  counter-propagating  modes 

u(z)  =  a(x)  -f  b(ar)  , 

where 


(44a) 
(44b) 

It  is  easy  to  confirm  that  these  modes  partition  both  the  power  flux  and  the  energy  density, 

P(a  +  b)  =  P(a)  +  P(b) 

E(a+b)  =  E(a)  +  E(b)  , 


a(x)  =  u_(a:o)  exp  [+i(x  -  i0)]  |~) 
b(x)  =  u+(xo)  exp  [-i(a:  -  x0)]  |+)  . 
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that  their  individual  energies  are  separately  conserved, 

P(a(x))  =  P(a(x0))  =  +c07|u_(a;o)|2 
P(b(a;))  =  P(  b(x0))  =  -co7|u+(x0)|2 

and  that  their  energy  densities  are  spatially  uniform, 

E(a(x))  =  E(  a(x0))  =  7|u_(x0)|2 
P(b(x))  =  £(b(x0))  =  7|u+(x0)|2  ■ 

This  is  d’Alembert’s  decomposition  again  —  an  exact  result  for  uniform  media,  but  only  approxi¬ 
mate  otherwise. 

4.  STATE-SPACE  TRANSFORMS 

Later  on,  the  convergence  of  Picard-series  solutions  will  be  improved  by  the  application  of  well 
chosen  transforms.  This  section  previews  the  properties  of  such  transforms  —  particularly  those  of 
the  wave-energy  conserving  “pseudo-unitary”  class. 

The  challenge  is  to  find  some  expression  for  the  solution  to  the  d’Alembert  initial-value  problem, 

dxK(x,xo)  =  G(x)K(x,xo)  ...  K(xo,xo)  =  l,  (45) 

that  is  better  behaved  than  Eq.  (37).  The  general  idea  is  to  do  this  by  transforming  to  still  other 
representations  where  the  generator  assumes  more  tractable  forms.  A  perfectly  diagonal  form  would 
seem  ideal,  but  it  is  not  generally  attainable.  Instead,  we  apply  a  series  of  transforms  that  improve 
the  diagonality  incrementally  with  each  step. 

The  first  step  is  along  the  general  lines  of  the  work  by  Keller  and  Keller  [6].  The  essence  of  it 
is  to  find  a  change  of  representation  that  converts  the  generator  G  into  the  form 

G  =  A  +  R  ,  (46) 


where 

•  A  is  diagonal  and  proportional  to  a  suitably  defined  quantity  v  «  n,  and 

•  R  is  off-diagonal  and  proportional  to  a  suitably  defined  small  quantity  q. 

Once  that  is  achieved,  one  can  transform  to  the  A-interaction  representation,  where  the  jfth  term 
in  the  propagator’s  Picard  series  will  contain  only  (p .  For  small  |£|,  this  approach  leads  to  a 
perturbative  formulation  for  the  field  in  terms  of  multiply  reflected  left-  and  right-going  modes. 
This  series  constitutes  an  algorithm  for  computing  the  field  at  limited  ranges,  and  its  leading 
term  provides  an  asymptotic  expression  for  the  field  at  long  range.  The  transform  that  splits  the 
generator  according  to  Eq.  (46)  is  a  simple  rotation  in  Pauli  space. 

In  fact,  the  second,  third  and  all  subsequent  steps  are  achieved  though  a  regular  succession  of 
such  rotations.  Given  certain  assumptions  about  the  range  dependence  of  the  medium,  each  one 
improves  both  the  perturbative  formulation  and  the  long-range  asymptotic  approximation.  The 
operators  that  produce  these  rotations,  however,  differ  in  two  respects  from  the  transform  U  in 
Eq.  (41),  namely,  they  are 
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•  ^-dependent,  and 

•  “pseudo-unitary”  (a  generalization  of  the  unitary  class). 

Ultimately  this  means  that  all  the  v  and  g  quantities  are  functions  of  x,  and  that  all  the  rotation 
angles  are  imaginary.  These  points  are  explored  in  the  remainder  of  this  section,  and  then  the 
rotations  themselves  are  obtained  explicitly  in  the  following  sections. 

Pseudo-Unitarity 

This  section  makes  some  analogies  to  quantum  mechanics,  and  thus  uses  the  Hamiltonian 
H  =  zG.  For  any  operator  V,  it  is  clear  from  Eq.  (42a)  that 

dx  (utVu)  =  uf  [i(HfV  -  VH)  +  d*v]  u  (47) 

represents  the  total  x-derivative  of  the  matrix  element  of  V  with  the  state  vector  u  —  the  rate 
of  change  due  to  both  the  dynamic  evolution  of  u  and  any  explicit  x  dependence  in  V.  Since  the 
Hamiltonian  is  “pseudo-Hermitian” ,  i.e. 

H*P  =  PH  . . .  pseudo-Hermitian  ,  (48) 

and  P  itself  is  independent  of  x,  any  matrix  element  of  P  must  also  be  x-independent.  That 
fact  itself  is  nothing  new  (it  was  implicit  in  Eq.  (19a));  the  novel  element  is  its  relation  to  the 
pseudo-Hermiticity  of  H.  Note  that  there  is  no  comparable  result  for  E.  Since 

E  =  -PH  ,  (49) 

that  operator,  although  pseudo-Hermitian,  is  x-dependent  wherever  H  is.  As  a  result,  P  is  a 
constant  of  the  motion,  but  E  is  not. 

The  dynamic  evolution  of  the  state  vector  can  be  summarized  in  terms  of  the  propagator  as 
u(£)  =  K(£,£)u(£).  From  the  identity 

uf(OPu(o  =  ut(o  [Kt(*,oPK(e,o]  u(o , 

together  with  the  initial  value  K(£,£)  —  1,  it  follows  that,  since  P  is  an  invariant  of  the  motion, 
the  propagator  is  “pseudo-unitary”,  i.e. 

P  =  Kt(£,£)PK(£,£)  ...  pseudo-unitary  (50) 

for  any  £,  £.  This  property  is  responsible  for  the  invariance  of  P  during  the  evolution.  By  definition, 
P  will  be  invariant  under  any  pseudo-unitary  transform  M,  whether,  like  K,  it  is  related  to  the 
dynamic  evolution  of  the  system  or  not.4 

Using  the  fact  that  P  is  Hermitian,  the  pseudo-Hermiticity  of  H  and  the  pseudo-unitarity  of  a 
transform  M  can  be  rewritten  in  various  equivalent  ways,  such  as 

(PH)*  =  PH 
(PM)*  =  PM-1  . 

4Strictly  speaking,  the  way  things  have  been  defined,  it  is  M-1,  not  M  itself,  that  is  directly  analogous  to  K. 
However,  pseudo-unitarity  for  one  implies  pseudo-unitarity  for  the  other. 
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These  results  have  better-known  analogs  in  quantum  mechanics,  where  the  operator  P2  =  1  is  a 
constant  of  the  motion  rather  than  P  itself.  Although  the  vibrating  string  problem  and  the  two- 
state  quantum  problem  are  very  similar  in  a  mathematical  sense,  it  would  be  a  mistake  to  press  the 
analogy  too  far,  since  the  two  state  vectors  have  entirely  different  physical  meanings  —  one  being 
a  classical  displacement  and  the  other  a  probability  amplitude. 

In  a  quantum  mechanical  setting,  the  physical  observables  are  represented  by  Hermitian  oper¬ 
ators.  Any  such  operator  V  has  real  Pauli  coefficients  Vo,  •  •  • ,  V3.  In  the  present  context,  V  may 
be  pseudo-Hermitian,  in  which  case  Vo  and  V3  are  real  while  Vi  and  V2  are  imaginary. 

State-Space  Transforms 


A  linear  transform  Mu  =  u  will  preserve  the  appearance  of  the  problem 


dx  u  =  Gu 

(51a) 

P(  u)  =  co7U*Pu 

(51b) 

E(u)  =  7U^Eu 

(51c) 

in  terms  of  the  transformed  operators 

G  =  A  +  R 

(52a) 

P  =  MfPM 

(52b) 

E  =  M*EM 

(52c) 

The  first  part  of  the  generator, 

A  =  m!gm 

(53a) 

is  a  similarity  transform;  the  second  part, 

R  =  — 

(53b) 

involves  the  explicit  x  dependence  of  M  (which  is  assumed  to  be  nonsingular  and  differentiable  so 
that  R  is  well-defined). 

To  conserve  energy,  M  will  be  restricted  to  pseudo-unitary  transforms.  When  that  is  done,  Eq. 
(49)  is  enough  to  determine  the  forms  of  the  transformed  energy  operators: 

E  =  -iPA  (54) 

and  P  =  P.  In  quantum  mechanical  terminology,  such  an  M  is  an  innocuous  “change  of  picture,” 
and  all  that  is  needed  to  characterize  its  effect  is  usable  expressions  for  A  and  R.  To  obtain  them, 
it  will  be  sufficiently  general  to  use 


M  =  exp(— iS)  (55) 

with  a  pseudo-Hermitian  S  to  guarantee  that  M  is  pseudo-unitary.  Thus, 


A  =  exp(+iS)G  exp(— iS) 
R=  —  exp(-fiS)cfa;exp(— iS)  . 


(56a) 

(56b) 
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Since  Eq.  (56a)  is  a  similarity  transform,  A  =  Aol  +  A-<x  is  the  result  of  a  simple  Pauli-space 
rotation  of  G,  as  discussed  in  Appendix  A.  Specifically, 

A0  =  G0  =  0  (57a) 

and  A  =  AA  where 

A  =  G  =  in,  (57b) 

with  A  obtained  by  rotating  g  about  the  s  direction  through  the  angle  <f>  =  -25.  But,  since  S  is 
pseudo-Hermitian  rather  than  Hermitian,  this  angle  need  not  be  real.  If  s  points  mainly  in  the  e3 
direction  (i.e.,  if  Js-e3 12  >  |s-ei|2  +  |s-e2|2),  then  <f>  is  &  real  angle  —  otherwise  it  is  imaginary. 

The  second  term  of  the  transformed  generator  has  the  form 

R  =  -  |e+t5°  [cos(5)  +  isin(5)s-<r]|  dx  je-iS°  [cos(5)  -  i  sin (5) ser]}  .  (58) 

With  the  help  of  s-s  =  0  and  Eq.  (A6),  this  reduces  to 

R  =  iSol  +  i  [5  s  +  cos(5) sin(5)s  —  sin2(5)(sxs)J  a  .  (59) 

(As  always,  the  dot  indicates  differentiation  by  the  scaled  coordinate  x.)  The  scalar  part  of  S 
affects  only  the_  scalar  term  Rq  =  iSq.  Since  So  has  no  impact  on  the  vector  part  of  R  and  had 
none  at  all  on  A,  it  is  entirely  irrelevant  to  the  problem  of  diagonalizing  A  +  R.  For  convenience, 
let  5o  =  0  so  that  only  the  vector  part 

R  =  —  -  [^>s  +  sirups  +  (1  —  cos<£)sxsj  (60) 

remains.  This  involves  rates  of  change  for  both  the  length  and  direction  of  5  (i.e.,  both  <j)  and  i). 
and  that  appears  to  be  about  all  that  can  be  said  about  it  in  general.  However,  if  the  direction  is 
fixed,  the  expression  simplifies  to 

R  =  ■  (61) 

One  case  is  useful  enough  to  deserve  special  attention,  namely  an  5  that  is  parallel  to  one  of  the 
Pauli  axes:  s-ej  =  ±1  for  j  =  1,  2,  or  3.  To  maintain  pseudo-unitarity,  (j)  must  be  real  when  j  —  3 
and  imaginary  when  j  =  1,2. 

5.  BREMMER  REPRESENTATION 

This  section  implements  a  pseudo-unitary  transform  to  the  “Bremmer"  representation,  where 
the  generator  of  the  wave  evolution  is  diagonal  to  order  e.  This  fact  improves  the  convergence 
of  the  Picard  series  for  the  solution  and  provides  the  impetus  for  additional  transforms  that  are 
developed  in  the  next  section  to  improve  it  still  further. 
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Change  of  Notation 

The  previous  section  dealt  with  state-space  transforms  in  a  fairly  general  sense:  u  =  Mu,  where 
u  was  a  d’Alembert  state  vector  and  u  was  the  state  vector  in  some  transformed  representation. 
Here  we  devise  a  particular  form  for  that.  Since  it  will  ultimately  be  the  first  in  a  series,  we  switch 
to  an  indexed  notation  at  this  point.  Quantities  in  the  d’Alembert  representation  together  with 
the  transform  that  produces  them  are  indicated  by  a  superscript,  e.g.,  G,  A,  R,  u  and  M, 
S  are  written  hereafter  as  G^°\  A^,  R(°\  and  M^,  S^.  The  transformed  or  “Bremmer” 
representation  is  labeled  with  a  ^  superscript,  so  G,  A,  R,  u  are  denoted  G^\  A^\  R^\ 

The  transform,  then,  is  written 


Rotation 


Making  A(1)  diagonal  is  tantamount  to  aligning  with  £3.  Since  A^1)  results  from  a  rotation 
of<?(°)  =  ee2— i(l+e)e3,  the  diagonalization  can  be  accomplished  by  using  e\  as  the  axis  and  rotating 
through  the  angle  tan-1(ie/(l  +  e)),  i.e.,  by  using  a  vector  with  components  =  0 

and 

=  —  «/>(0)/2  ,  (62) 

where 

tanhV>^  =  —  •  (63) 

Thus 

S(°)  =  -|v-(0Vi  ,  (64) 

and  the  transform  =  exp  {— ;s(°)}  is 


M(0) 


The  transformed  operators  are 


where 


=  cosh  /2j  1  —  sinh  / 2^  <T\ 

1  +  ni  t  1  —  n 
~  2^H1  1  2y/n <Tl 

(65) 

H- ‘ 

II 

1 

c«». 

t— » 

CO 

(66a) 

R(J)  =  gWo-i 

(66b) 

CO 

b 

1 

II 

*-4 

(66c) 

E(1)  _  , 

(66d) 

jn  =  n 

(67) 

£>(1)  =  =  d*  0°g  • 

(68) 

Note  that  is  precisely  the  condition  anticipated  by  Eq.  (31)  at  the  end  of  section  2. 
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Using  a  pseudo-unitary  transform  automatically  produces  an  off-diagonal  R/1).  The  converse 
is  also  true.  Keller  and  Keller  [6]  began  by  requiring  R/1)  to  be  off-diagonal  and  arrived  at  a 
pseudo-unitary  transformation  (the  same  one,  to  within  an  arbitrary  constant).  Pseudo-unitarity 
distinguishes  an  essentially  unique  “splitting”  among  all  the  possible  transforms  [7].  It  is  appro¬ 
priate  to  acknowledge  here  that,  despite  a  somewhat  altered  appearance  in  the  present  context, 
the  transform  developed  in  this  section  is  a  type  of  Foldy-Wouthuysen  transform  [8].  This  sort  of 
transform  was  introduced  in  the  early  1950s  [9]  to  produce  asymptotic  mode  decoupling  in  a  very 
different  physical  context  —  relativistic  quantum  mechanics. 

Mode  Separation 

This  is  a  convenient  point  to  pause  and  consider  the  case  of  an  environmental  region  where  the 
wavenumber  k  is  a  constant  —  possibly  different  from  the  reference  value  fco-  Since  e  is  constant 
in  such  a  region,  G^0)  is  independent  of  x  and  can  thus  be  diagonalized  by  an  ^-independent 
This  means  that  R^1)  vanishes  so  that  Eq.  (64)  diagonalizes  the  whole  generator  =  —  zi/1)  <73, 
and  that  =  n  is  independent  of  x.  As  in  section  3,  the  generator’s  eigenvectors  are  |c);  here, 


however  the  corresponding  eigenvalues  are  —  Again  the  state  vector  separates  exactly  into  a 

sum  of  counter-propagating  modes 

u^(x)  =  a^(x)  +  b^(x)  , 

(69) 

where 

a^(x)  =  it^(xo)  exp  [+zV^  x  (x  —  xo) 

I-) 

(70a) 

b^(x)  =  vf+\x, o)  exp  x  (x  —  xq) 

l+). 

(70b) 

It  is  easy  to  confirm  that  these  modes  also  partition  both  the  power  flux  and  the  energy  density, 

P(  a(1)  +  b(1))  =  P(a(1>)  +  P(  b(1>) 

E( a(1)  +  b(1))  =  E( a(1>)  +  E{ b(1))  , 


that  their  individual  energies  are  separately  conserved, 


P(a(1)(x))  =  P(a(1)(x0))  =  +co7|u^(zo)|2 
P(b(1\x))  =  P(b(1)(x0))  =  -cot|u+Vo)|2  , 


and  their  energy  densities  are  spatially  uniform, 

E(  a^(x))  =  P(a^(xo))  =  7^^|u^(xo)|2 
P(b(1)(x))  =  E( b(1)(x0))  =  7t'(1)ju+)(xo)|2  • 


This  is  the  d’Alembert  decomposition  of  section  3  again,  this  time  for  if  n  ^  1. 
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Bremmer  Series 

In  the  Bremmer  representation,  the  initial-value  problem  is 

dxU^(x)  =  jyl^^x)  +  ItP^x)]  u^(x)  ...  u^^(xo)  =  [m^°^(xo)]  u^°^(xo)  . 

Since  il^(  x),  the  “large”  part  of  the  generator,  is  diagonal,  it  can  be  used  to  construct 
action”  transform  [10], 

L(x,x0)  =  exp  d£A(1)(£)^  = 

in  which 

¥>(1)(x,x o)  =  /  dfi/(1)(0- 

JXQ 

This  L  converts  the  state  vector  to  the  interaction  representation  via 

u(1)(x)  =  L(x, x0)q(x)  , 

and  is  itself  the  propagator  for  an  artificial  problem, 

dxL(x, xo)  =  yl(1)(x)L(x, xo)  ...  L(xo,  xo)  =  1 
in  which  the  perturbation  g ^  is  “turned  off”  [2]  so  that  |+)  and  |— )  are  decoupled. 

In  this  interaction  representation,  the  initial-value  problem  for  state  vectors  is 
dxq(x)  —  V(x,x0)q(x)  ...  q(x0)  =  u(1)(xo) 
in  terms  of  the  generator,5 

V(x,xo)  =  L_1(x, xo)R/^(x)L(x,  xo)  . 

The  propagator,  Z(x,xo),  evolves  the  state  vector  from  its  initial  value  in  the  usual  way, 

q(x)  =  Z(x,  x0)q(x0)  , 

and  thus  satisfies  the  initial-value  differential  problem, 

dzZ(x,x0)  =  V(x,x0)Z(x,x0)  ...  Z(x0,x0)  =  l. 

As  before,  the  equivalent  integral  equation  is 

Z(x,  x0)  =  1  4-  [  d£  V (£,  xo)Z(£,  x0)  , 

Jx  o 

and  its  Picard  series  solution  (often  called  a  Born  series  in  this  context  [10])  is 
5In  quantum-mechanical  terms,  the  interaction  Hamiltonian  would  be  iV. 


exp 

0 


exp  f+t^^aqxo)) 


(71) 
“inter- 

(72) 

(73) 

(74) 

(75) 

(76) 

(77) 

(78) 

(79) 

(80) 
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00 

Z(x,x0)  =  j(x,x0)  , 

i= o 


(81a) 


where 


Zq(x,xo)  =  1 


(81b) 


and 

Z,j+i(x,x0)=  f  d£  V(^,x0)  Zj(f,x0)  • 

Jx  0 


(81c) 


Equations  (79)— (81)  are  identical  in  form  to  Eqs.  (35)-(37).  Here,  however,  V  a  gW; 


y(t,xo)  =  e(1)(0 


0  exp  [+2ty>W(£,xo)| 

exp  [-2t>(1)(^,x0)j  0 


(82a) 


Thus  Z j  involves  only  the  jth  power  of  so  the  convergence  problem  that  was  caused  by  the 
presence  of  terms  containing  both  e°  and  e1  in  G  (the  generator  in  Eq.  (39a))  is  avoided.  For 
|e(1)|  <C  1,  the  series  should  converge  rapidly. 


Furthermore,  V  is  off-diagonal.  Consequently  the  Z j  terms  are  alternately  diagonal  and  off- 
diagonal.  Z12,  for  instance,  is  just  the  double  integral  of 

V(£,xo)V(C,x0)  =  6{1)(0e(1)(0 

As  a  result,  when  Eq.  (81a)  is  used  with  Eq.  (78)  with  the  state  vector  expanded  in  the  unperturbed 
eigenstates,  q(x)  =  X)«=±  9<r(x)K)!  one  finds 


exp  [+2^(1)(^,C)J  0 

0  exp  f-2^1)(^,C)l 


(82b) 


9<(x)  =  9<;(xo)  +  J  dx\  g(1)(xi)exp  [+2<r^(1)(xl, xo)]  g_?(x0) 

+  J  dx2  J  dxi  g^l\x2)g{1\xi)exp  [+2^V(1)(x2,  xi)]  gt(x0) 


+ 


(83) 


for  <;  =  ±.  This  is  essentially  the  Bremmer  series  [11]  —  an  expression  whose  jth  term  represents  a 
wave  that  been  reflected  j  times  on  its  way  from  xq  to  x,  with  acting  as  a  distributed  reflec¬ 

tion  coefficient  for  the  medium.  Clearly,  if  vanishes  at  xo,  the  odd-numbered  terms  all  vanish. 
Then  the  state  vector  reduces  to  <fr(x)|<r),  a  mode  that  propagates  purely  in  one  direction,  but  has 
an  amplitude  gc(x)  composed  of  contributions  that  have  undergone  an  even  number  (0, 2,4,  •••) 
of  distributed  reflections.  The  j  >  2  terms  were  left  implicit  in  Eq.  (83).  More  of  them  could 
have  been  included  without  much  trouble,  but  there  would  be  little  reason  for  it  in  this  context 
because  our  interest  is  in  long-range  propagation.  For  computational  purposes,  the  infinite  series 
must  be  truncated  at  some  order,  and  this  implies  a  maximum  x  beyond  which  the  truncated  series 
no  longer  faithfully  represents  g?(x).  For  long-range  use,  a  large-a;  asymptotic  expansion  is  needed 
instead.  The  first  term  in  the  asymptotic  expansion  of  the  Bremmer  series  for  small  is  just  the 
leading  term  in  the  series  [6], 
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u{^(x)  ~  exp  q<(xo)  .  (84) 

This  expression  is  an  asymptotic  approximation  for  one-way  propagation  that  excludes  backscatter 
in  all  forms  (all  the  j  >  0  terms). 

Note  that  this  would  all  have  been  equally  simple  had  R(°)  involved  <7 2  instead  of  <T\,  i.e.,  if 
had  turned  out  to  be  orthogonal  to  e 2  rather  than  to  e\.  In  that  case,  a  diagonal  A/1)  could 
have  been  produced  by  simply  rotating  around  that  axis  instead. 

6.  HIGHER-ORDER  REPRESENTATIONS 

Building  on  the  results  of  the  preceding  section,  this  section  introduces  a  hierarchy  of 
trans-Bremmer  rotations.  In  these  the  wave-evolution  generator  is  diagonal,  and  thus  the  for¬ 
ward/backward  modes  are  decoupled,  to  progressively  higher  orders. 

Chain  of  Rotations 

Pseudo-unitary  rotations  similar  to  the  one  in  the  preceding  section  will  be  used  to  generate 
further  representations,  which  will  be  denoted  by  superscripts  (m),  with  m  >  1.  Let  ej ,  e#  be  the 
first  two  Pauli  unit  vectors  ei,e2  in  either  order  and  suppose  that  is  orthogonal  to  ej  (i.e., 
—  0).  Then  choosing  to  lie  along  the  j  axis  (|s(m)-ej|  =  1)  allows  the  first  term  on  the 
right-hand  side  of 

d(m+1)  =  (cosht/^G^  +  isinhV’(m)G(m)xs(m))  +  \s^dx^  (85) 

(see  Eq.  (All))  to  be  aligned  with  e3  by  the  proper  choice  of  This  is  still  only  a  partial 

diagonalization  because  the  second  term  remains  orthogonal  to  £3;  however,  if  |£'(TO+1)|/|£'(m+1)|  < 

|G^m^|/|G3  |,  it  is  a  step  in  the  right  direction.  Furthermore,  since  =  0,  the  process  can 

easily  be  repeated,  this  time  by  a  rotation  about  the  l  axis.  In  fact,  it  can  be  iterated  indefinitely 
in  an  alternating  series  of  rotations  about  the  first  and  second  Pauli-space  axes. 


The  process  begins  with  the  Pauli  vector  part  of  the  d’Alembert  generator  rewritten  as 


<?(°)  =  -i^(0)e3  +  0(O)  e3xs(0) 

(86a) 

in  terms  of 

S<°>  =  e\ 

(86b) 

«/<0>  =  l  +  e 

(86c) 

^°)  =  e. 

(86d) 

The  axis  and  the  rotation  angle 

^°>  =  tanh-1  (<?((V0)) 

(86e) 

determine  the  transform  that  produces  G^1). 
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The  jth  iterate  has  the  same  form, 


=  —iu^e 3  -f  ez  , 

(87a) 

with 

(87b) 

j/W)  =  cosh 

(87c) 

eU)  =  \dA{i-l) . 

(87d) 

The  axis  and  rotation  angle 

=  tanh-1 

(87e) 

determine  the  transform  that  produces  G^+1\ 

This  procedure  produces  a  chain  of  generators 

q(0)  _►  qC1)  _+  q(2)  — >  G^  —>•••—» 

(88) 

linked  by  a  series  of  pseudo-unitary  transforms,  where 

=  cosh  (pip^/ 2)  1  —  sinh  (if>^/2 )  s^-a  . 

(89) 

Only  one  assumption  is  implicit  in  the  construction  of  link  m  at  the  end:  that  exists  for  the 

environment  in  question.  As  m  is  increased,  lengthening  the  chain,  p(m)  always  remains  identical 

to  the  original  p(°\  and  the  system  evolves  by  the  transformed  equation 

dx u(m)  =  G(TO)u(m)  . 

(90) 

As  m  is  increased  by  1,  s simply  rotates  by  —  ir/2  about  £3.  (It  is  periodic: 

s(m+4)  =  s("0.)  In 

addition,  provided  the  ratio  |i/TO)/f>(ra)|  diminishes  with  increasing  m,  the  generator  becomes  more 

and  more  diagonal.  Together  these  mean  that  successive  G'm'  vectors  spiral  in 
with  —ez  as  sketched  in  Fig.  (1). 

toward  alignment 

As  in  the  m  —  1  case  of  section  5,  the  transformed  generator  for  m  >  1  is 
A +  R(m\  where 

split  into  G(m)  = 

A(rn)  = 

(91a) 

R(m)  =  _ 

(91b) 

And,  as  always, 

C*5 

b 

1 

II 

? 

(91c) 

Thanks  to  the  nature  of  the  transform  that  produced  them,  all  three  operators  have  conveniently 
simple  forms:  p(m)  is  invariant,  oc  is  diagonal,  and  oc  involves  only  one  of  the 
two  off-diagonal  Pauli  operators  0’\,0‘2-  It  would  be  too  much  to  hope  that  the  simplicity  of  Eq. 
(66d)  would  persist  in  a  form  like  =  z/(m)l,  and  indeed  it  does  not.  But  there  is  no  physical 
reason  why  it  should. 
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Fig.  1  -  Pauli-space  illustration  of  G(m*  for  m  =  0  •  •  •  4.  Since  G(m*  is  imaginary,  the  figure  is  only  qualitative.  The 

mth  gray  arrow  has  components  (Gjm\  G^,  QG(m)). 


Mode  Separation 

Suppose,  for  the  moment,  that  the  environment  has  an  interval  I  where  g^m\x)  =  0.  Then  the 
generator  is  diagonal  in  that  interval 

G^m^(x)  =  —  ii/(m\x)(T3  ...  x  &  I 


and  so  is  the  propagator 

K(m)(x,x0)  =  exp  J  d<i/(m)(C)j  ...  x,x0El, 

so  that  the  state  vector  separates  exactly  into  a  sum  of  counter-propagating  modes 

u(m)(x)  =  a(m>(x)  +  b(mHx)  , 

where 


sSm\x)  =  u^™\xo)exp^+i  J  ()  |— ) 

b(m\x)  =  u+'\x0)exp^-i  J  d(u^m\()  |+)  . 


(92) 

(93a) 

(93b) 
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It  is  easy  to  confirm  that  the  power  flux  is  partitioned  between  these  modes 

P(  a(TO)  +  b(m))  =  P(  a(ra))  +  P(b(m)) 

and  separately  conserved 

P(a(m)(x))  =  P(a(m)(x0))  =  +co7|«(-^‘)(a;o)|2 
P(b(m)(*))  =  P(b(mW)  =  -co7|«(+ra)(xo)|2  . 


This  result  is  exact  but,  of  course,  it  is  not  general  because  the  condition  Q^m\x)  —  0  limits 
its  validity  to  a  particular  type  of  environment.  The  m  =  1  case  was  encountered  in  section 
5,  where  ^(x)  =  0  was  seen  to  imply  a  uniform  medium  having  n(x)  =  n(xo),  with  standard 
d’Alembert  mode  separation  prevailing  everywhere.  The  limitation  is  even  greater  for  m  >  1. 
For  example,  g^\x)  =  0  presupposes  a  medium  where  n~l(x)  -  n(xo)-1  =  (x  —  xo)g  for  some 
constant  g.  This,  in  turn,  means  that  0  <  n(x)  <  oo  is  possible  only  on  the  semi-infinite  interval 
I  =  {x|(x  —  xo)gn(xo)  >  0}.  This  is  just  another  reminder  that,  except  for  tailor-made  special 
cases,  mode  separation  is  an  approximate  result,  not  an  exact  one. 

Higher-Order  Bremmer  Series 

In  parallel  to  section  5,  the  initial- value  problem  for  u  can  be  solved  via  the  A interaction 
representation.  The  eigenvectors  of  A are  still  |c),  and  the  counterpart  of  Eq.  (83)  for  the 
components  of  the  interaction-representation  state  vector  is  again  a  Bremmer  series, 

q<(x)  =  Qc(xo)  +  J  dxi  £(m)(xi)exp  [+2^(rn)(a;l> ^o)]  9-c(*o) 

+  [  dx 2  f  dxi  Q('rn'> (X2) Q^ixi)  exp  \+2<;ig^m\x2,  aq)]  qq(xo)  H -  (94) 

Jxo  Jxo  1  1 


The  first  term  in  the  asymptotic  expansion  of  this  mth-order  Bremmer  series  for  small  is  again 
just  the  leading  term  in  the  series, 

u^m)(x)  ~  exp  [-<;zV(m)(:E,xo)]  g?(x0)  ...  C  =  ±  .  (95) 

This  is  a  higher-order  asymptotic  approximation  for  one-way  propagation. 

7.  EXPANSION  IN  e  DERIVATIVES 

This  section  introduces  a  smallness  criterion  for  spatial  derivatives  of  e  that  renders  the  mth- 
order  Bremmer  representation’s  generator  diagonal  to  order  em,  thus  eliminating  backscatter  to 
that  order.  Further,  it  obtains  the  resulting  x-dependent  amplitude  and  phase  of  the  field  and 
interprets  each  of  these  as  having  a  contribution  that  accumulates  along  the  propagation  path  and 
another  contribution  from  the  path’s  endpoints. 

If  the  environmental  variations  are  adiabatic  in  the  sense  that  x  ~  e-1,  then  one  should  expect 

that 


d™e  ~  em+l 


(96) 
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Since  tanh  =  e/(l  +  e)  and  tanh  if>^  =  dxe  /  2n,  one  would  have 

tanh  rv  £ 

tanh  (V  £^ 


and,  indeed,  for  any  order  m, 


tanh^(m)  ~  em+1  .  (97) 

The  symbolic  mathematics  system  Maple®  has  been  used  to  confirm  through  m  =  5  that,  assuming 
Eq.  (96), 

tanh^m)  =  2 +  0(em+2)  , 

which  provides  a  direct  verification  of  Eq.  (97).  This  means  that,  as  m  increases,  g(m)  grows  steadily 
smaller  in  relation  to  v in  the  sense  that  ~  em+1.  Thus,  under  these  assumptions, 

G^m)  is  diagonal  to  order  em. 


The  notation  has  been  used  generally  for  any  quantity  (  in  the  mth  representation.  When 
the  e  series  for  is  truncated  at  the  mth  order,  the  result  will  be  denoted  £H.  For  the  generator 
G(m),  this  is 


GH  =  -jvWff 3  ,  (98) 

where  i/M  consists  of  all  the  terms  of  up  through  order  em,  i.e.,  =  i/H  +  0(em+1). 

Since  Eq.  (98)  is  diagonal,  the  evolution  that  it  generates  can  be  found  by  a  trivial  integration, 
provided  that  the  function  v^m\x)  is  known.  As  m  grows  larger,  the  calculation  of  i/H  becomes  an 
increasingly  intricate  task,  but  one  that  reduces  to  a  fixed  pattern  of  routine  operations  —  an  ideal 
job  for  symbolic  computation  software.  Maple  has  been  used  to  generate  the  result  out  to  m  =  6. 
(Farther,  actually,  but  this  has  to  be  stopped  somewhere.)  The  result  can  be  expressed  as 


1 

+r 


3 


-  (r4 + r2) +  + H  -  + lA2 + k2) +  ° «7>  ■ 


(99) 


with  i/(ml  found  by  simply  truncating  v  at  the  mth  order.  Maple  has  also  been  used  to  verify  that 

eH  =  0  =  ^H. 


Cumulative  Contribution 

With  1  (x)  in  hand,  it  is  a  trivial  matter  to  solve  the  problem  approximated  at  order  m, 

d*u(m)  =  GHu("0  .  (100) 

The  state  vector  has  the  separated  form 


where  the  modes  are 


u(m)  _  a(m)  +  ^(m)  , 
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a(m)(x)  =  u+^xo)  exp  (+i^mkx,xo))  |— ) 

(101a) 

b(m)(x)  =  u^\x o)  exp  iy?lml(x,x0)^  |+) 

(101b) 

in  terms  of  the  phase  function 

<p(x,x o)  =  [  d£z/(£)  . 

Jx  0 

(102) 

(<pM  is  formed  in  the  usual  way,  by  using  z^m'  as  the  integrand.)  A  few  points  should  be  noted 
here: 

1.  The  important  quantity  here  is  this  phase,  specifically  the  deviation  from  y?[°l(x,xo)  =  x  —  xq 
(its  value  in  a  uniform  medium  with  n  =  1). 

2.  The  even-order  contributions  to  are  all  negative.  For  x  —  xq  >  0  (<  0),  this  means 

that  all  the  even-order  contributions  to  serve  to  retard  (advance)  the  phase  in  a  strictly 
cumulative  way. 

3.  The  odd-order  contributions  to  do  not  have  definite  signs  (thanks  to  their  odd  powers  of 
e)  so  they  are  likely  to  do  quite  a  lot  of  zero-crossing  when  |x  —  xo|  is  large  enough.  In  that 
event,  their  contributions  to  should  oscillate  about  zero  rather  than  accumulating. 

4.  Comparison  of  Eq.  (99)  to  the  Taylor  series  for  the  refractive  index, 

,  1  2  1  3 

n  =  l  +  e--e  +  -e3 

5  4  7  c  21  c  .  y. 

-8e  +8e  -ie'  +0(e  ) 

shows  that  v  begins  to  deviate  from  n  at  order  e4.  Thus  the  WKB  phase  approximation  —  the 
consequence  of  invoking  the  estimate  u  «  n  in  the  phase  integrand  —  begins  to  depart  from 
the  above  results  at  fourth  order. 

Endpoint  Contribution 

The  procedure  detailed  above  yields  u.(m\  the  state  vector  in  the  mth  transformed  representa¬ 
tion.  That  still  has  to  be  transformed  back  to  the  d’Alembert  (0-th)  representation  via 

u(0)  =  W(m)u(m)  1  (103) 

where 

W(m)  =M(0)M(1)---M(m_1)  .  (104) 

So  W<m)  needs  to  be  calculated.  Of  course,  since  it  is  going  to  be  applied  to  the  mth  order 

estimated  state  vector  that  results  from  approximating  G^7"^  by  G  t,nl ,  a  suitable  approximation 
for  will  do.  But  what  is  “suitable”?  A  knee-jerk  option  might  be  Wlml,  the  transform 

truncated  at  the  same  order  as  the  generator.  On  consideration,  however,  that  appears  to  be 
overdoing  it.  The  phase  deviation  ip ^  embodies  the  effects  of  environmental  nonuniformities 

accumulated  throughout  the  propagation  from  £  =  xo  to  £  =  x.  But  the  back-transform  is 

not  cumulative;  it  depends  only  on  environmental  conditions  at  the  final  point  £  =  x.  It  might  be 
sufficient  to  approximate  it  by  for  some  l  <  m  .  We  proceed  for  the  moment  with  WM. 


Systematic  Splitting  of  Wavefields 


25 


The  unitary  transform  in  Eq.  (41)  needs  to  be  applied  to  return  to  the  initial  representation 
(where  the  first  component  of  the  state  vector  was  simply  the  string  displacement  and  the  second 
was  the  slope  of  its  tangent) 


a  =  Ua(0)  =  UWMa(<)  (105a) 

and 

b  =  Ub(0)  =  UW^b(€)  .  (105b) 


In  view  of  Eqs.  (93)  and  (101),  only  the  matrix  elements  (+|UWl£l |+)  and  (+ |UW^|— )  are  needed 
to  obtain  the  displacement  fields  for  the  left-  and  right-going  d’Alembert  modes,  respectively. 
Another  resort  to  Maple  yields 

(+|UW^|+)  =  (106a) 

(+|UWM|-)  =  ,  (106b) 

V2 

where  the  phase  is 

t?  =  —  fee  +  (^e2e  4-  pje)  +  0(c5)  ,  (107) 


and  the  amplitude  term  is 


Q  = 


(108) 


The  endpoint  phase  vanishes  to  first  order,  i?!1'  =  0.  When  the  endpoint  x  lies  in  a  uniform  region 
(where  n  is  constant,  but  not  necessarily  1),  the  phase  vanishes  entirely,  and  if  n  =  1  the  amplitude 
term  also  reduces  to  unity.  Comparison  of  Eq.  (108)  to  the  power  series 


n",/2=1-r+r2 

15  q  195  4  c. 

-16e  +S£  +0(t) 


(109) 


shows  that  the  standard  WKB  amplitude  expression  Q  «  l/\/n  is  valid  only  through  order  e2.  It 
might  be  guessed  that  Q  «  1  /y/u  would  be  an  improvement  on  that,  but  the  e-series 


shows  that  it,  too,  is  valid  only  through  order  e2. 
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Conclusion 

The  final  result  for  the  string  displacement  field  is 

u(x)  =  (u(Hm)(a:o)exp|+i[^[ml(x,xo)  -  tf[£1(x)]}  +  u^T^o)  exp  (x ,  x0)  -t?M(x)]})  , 

(110) 

where  m  is  the  “cumulative”  order  and  l  <  m  is  the  “endpoint”  order. 

8.  EXAMPLE 

This  section  introduces  a  specific  weakly  nonuniform  medium  and  presents  a  numerical  example 
of  4th-order,  non-WKB  phase  accumulation. 


Consider  a  class  of  environments  characterized  by  a  periodic  inhomogeneity  of  the  form 

e(x)  =  Ag  sin(agx)  +  Bg  sin(bgx)  .  (Ill) 

As  an  example,  take  the  case  where  (A,  B,a,b,  g)  =  (1.0, 0.4, 1.0, 0.8, 0.1).  Figure  2  illustrates  the 
first-  through  fourth-order  contributions  to  Eq.  (99),  namely, 


first-order  index  contribution - he 

1  2 

second-order  index  contribution - -e 

third-order  index  contribution - h^e3 

z 

5  a  1  2 

fourth-order  index  contribution  ■■■—-€  — -e 

WKB  non- 
WKB 


plotted  from  top  to  bottom  in  the  figure.  As  should  be  expected,  this  e(x)  exhibits  spatial  inter¬ 
ference  in  the  form  of  a  pattern  of  beats  in  the  top  plot  with  a  spacing 

L  «  27r/[(a  —  b)g)  =  1007T  , 

and  the  second  through  fourth  orders  follow  suit.  The  amplitude  of  the  beats  in  the  mth  order  is 
gm  =  10_m,  with  the  odd  orders  being  about  equally  positive  and  negative,  and  the  even  orders 
purely  negative.  In  the  bottom  plot,  the  small  non-WKB  part  of  the  fourth-order  contribution  is 
shown  along  with  the  total. 


Figure  3  shows  the  phase  contributions  according  to  Eq.  (102).  Each  plot  results  from  inte¬ 
grating  the  corresponding  plot  in  Fig.  2  from  xq  =  0  to  x.  Their  behavior  is  as  anticipated:  the 
odd  orders  oscillate  around  small  ‘dc’  values;  the  even  orders  decrease  monotonically,  contribut¬ 
ing  steadily  accumulating  phase  retardations.  In  the  even  case,  the  phase  retardation  rate  drops 
sharply  with  increasing  order.  In  the  bottom  plot,  the  WKB  and  non-WKB  parts  of  the  fourth- 
order  contribution  are  shown  in  addition  to  the  total.  In  that  plot,  the  rate  of  non-WKB  phase 
retardation  is  about  15%  of  the  WKB  rate,  which  is  itself  only  about  1.5%  of  the  second-order  rate. 
Clearly,  fourth-order  effects  could  be  prominent  only  at  long  ranges.  In  fact,  it  is  only  when  ranges 
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-0.01  -> 
o 


Fig.  3  -  Phase  contributions  to  Eq.  (102)  for  the  example 
are  computed  using  20-digit  precision.  (The  tiny  ripples  ii 
the  same  when  computed  with  10-digit  precision.) 
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of  x  «  224, 000  ~  714 L  are  reached  that  the  fourth-order  non-WKB  phase  contribution  reaches 
— tt/2.  At  that  range,  the  WKB  contribution  to  the  fourth-order  phase  is  — 37r,  while  the  second- 
order  contribution  is  approximately  —  2147T  =  107  *  (— 2ir).  The  example  serves  to  underscore  the 
following  points: 

•  For  situations  that  demand  the  relative  phase  of  the  signal  at  points  only  a  few  dozen  wave¬ 
lengths  apart  (Ax  «  10  to  100),  fourth-order  phase  corrections  are  quite  small.  In  practice,  a 
first-order  estimate  will  usually  be  accurate  enough. 

•  When  the  demand  is  for  accurate  phase  values  at  very  long  ranges  (e.g.,  in  some  long-range 
interferometry  and  long  time-of-flight  applications),  fourth-order  phase  estimation,  including 
the  non-WKB  contribution,  may  be  indicated. 

9.  DISCUSSION 

For  the  one-dimensional  case  of  interest  here,  it  has  been  shown  that  a  properly  chosen  series 
of  pseudo-unitary  transforms  converts  the  equation  for  the  dynamic  evolution  of  the  wave  field 
—  the  Helmholtz  equation  in  this  case  —  to  a  fair  more  tractable  (in  fact,  diagonal)  approximate 
form.  This  can  be  done  to  any  desired  order  m  in  the  small  environmental  inhomogeneity  e(x). 
The  resulting  equation  can  be  integrated  immediately  to  provide  mth-order  expressions  for  the 
amplitude  and  phase  of  the  wave  field  at  long  range.  An  inherent  part  of  this  construct  is  that 
backscatter  is  neglected  at  all  orders  so  that  the  field  consists  of  a  pair  of  asymptotically  decoupled 
counter-propagating  modes. 

The  approach  taken  in  sections  6  and  7  is  not  without  historical  precedent.  Most  notable  is 
the  1950  application  by  Foldy  and  Wouthuysen  [9]  of  such  methods  in  quantum  physics  to  obtain 
the  Schrodinger  equation  from  its  relativistic  precursor,  the  Klein-Gordon  equation.  More  recently, 
Wurmser  et  al.  [12]  brought  these  techniques  to  bear  on  a  case  of  classical  wave  motion  in  which  the 
medium  is  also  allowed  to  vary  in  directions  transverse  to  the  direction  of  propagation  (e.g.,  along 
y  as  well  as  x).  The  resulting  two-dimensional  Helmholtz  equation  (Eq.  (32a)  with  — ♦  d%.  +  dy) 
leads  to  a  d’Alembert  representation  of  the  dynamics  that  is  similar  to  Eq.  (42a)  except  that  e  is 
a  differential  operator  in  y.  The  upshot  is  a  pair  of  asymptotically  decoupled  parabolic  equations 
for  propagation  in  the  ±x  directions.  These  involve  fourth-order  non-WKB  “corrections”  to  the 
refractive  index  like  those  seen  here.  They  also  contain  novel  features  related  to  the  additional 
transverse  degree  of  freedom,  notably  a  classical  analog  of  the  quantum  mechanical  phenomenon 
of  Zitterbewegung  [8]. 

Bremmer  originally  developed  the  representation  of  section  5,  Eq.  (83)  in  particular,  by  different 
means  [11].  He  first  “stratified”  the  medium,  approximating  the  refractive  index  by  a  piecewise 
constant  stairstep  function;  then  he  accounted  for  all  the  multiple  reflections  at  each  discontinuity; 
and  finally,  he  took  the  limit  of  infinitely  many  steps  of  vanishing  height.  This  “infinitesimal” 
approach  is  perfectly  valid  and  can  lead  to  physical  insights  [13].  It  has  been  avoided  here  only 
because  it  is  difficult  to  generalize. 
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Appendix  A 

PAULI  OPERATORS 


Pauli  matrices  are  used  to  facilitate  the  analysis  throughout  the  body  of  this  work.  This  section 
reviews  their  properties  —  especially  their  role  as  infinitesimal  generators  of  rotation  transforms  in 
state  space. 

The  Pauli  spin  matrices6 


(T<i  =  i 


(Al) 


have  seen  long  service  as  a  labor-saving  device  in  quantum  mechanics  and  will  prove  useful  here 
also.  Although  certainly  not  essential,  they  do  have  some  attributes  that  reduce  the  work  involved 
in  doing  linear  algebra  on  complex  two-dimensional  state  spaces,  e.g., 

trace  <Jj  —  0  =  1 

det<Tj  =  -1  <rj  =  <Tj. 

In  addition,  the  eigen-basis  of  <r3: 

<T3|<;)  =  <?k)  ...  C  =  ± 

where 


is  convenient  for  representing  complex  vectors 

q  =  H  • 

<=± 

The  main  advantage  of  Pauli  operators,  however,  relates  to  the  representation  of  linear  transforms. 

Pauli  Space 

The  three  Pauli  operators  are  complete  in  the  sense  that,  together  with  the  unit  matrix,  they 
form  a  basis  for  complex  2-by-2  matrix  operators.  Any  such  operator  can  be  represented  as 

V  =  F01  +  Vi*i  +  V2tr2  +  F3<r3  (A2a) 

8No  distinction  is  observed  between  the  operators  themselves  and  their  matrix  representation. 
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in  terms  of  the  four  coefficients 

Vo  =  trace  (V)/2  (A2b) 

Vj  =  trace(V<7j  +  c r,- V)/4  ...  j  =  1,2,3,  (A2c) 

all  of  which  have  real  values  if  V  is  Hermitian.  In  brief, 

V  =  VQ1  +  v  a  ,  (A3) 

where 

(T  =  (T (T 2&2  "t"  O’ 3^3  (A4a) 

V  =  Fiei  4-  V2S2  +  V^e3  (A4b) 


are  vectors  in  an  abstract  three-dimensional  ‘Pauli  space’.  The  original  operator  is  equivalent  to  a 
scalar  part  Vo  and  a  Pauli- vector  part  V.  A  magnitude  and  direction  can  be  defined  for  this  vector 
part  through  V  =  (F-F)1/2  and  v  =  F/F,  although  these  need  not  generally  be  real  valued. 

Products  and  Exponentials 

Products  and  exponentials  of  Pauli  operators  can  be  evaluated  by  using  the  fundamental  com¬ 
bination  rule 


Gjt Tfc  —  'i£jkt&l  "V  bjk  1 
and  its  various  spinoffs  [Al]  such  as 

(a-<r)(b-(r)  =  (a-b)l  +  i(axb)-cr 


and 


exp (dtziSs-a)  =  cos(S)l  ±  isin(5)s-<r  . 


(A5) 


(A6) 


(A7) 


Similarity  Transforms  as  Pauli-Space  Rotations 

For  any  Hermitian  operator  S,  the  result  of  using  the  unitary  operator  exp(— zS)  to  perform  a 
similarity  transformation 


V1  =  exp(+zS)V  exp(— z’S)  (A8) 

on  the  V  of  Eq.  (A3)  is  necessarily  an  operator  of  the  same  form, 

V'  =  F0'l  +  F'  <x  .  (A9) 

Clearly,  the  scalar  parts  of  S  and  V  play  only  a  trivial  role  in  this  (So  is  irrelevant  and  Fo  is 
invariant),  and  the  vector  parts  are  related  through 

F'-<r  =  exp(-t-z5-<7)F.<7  exp(— iS-(r)  .  (A10) 

Expansion  of  the  exponentials  using  Eq.  (A7),  followed  by  two  applications  of  Eq.  (A6)  results  in 


Systematic  Splitting  of  Wavefields 


33 


V'  =  cos(2S)  [V  -  (F-s)sJ  +  sin(2S)Yxl  +  (F-s)s  .  (All) 

To  clarify  the  form  of  this,  let  S  =  -</>/2  and  write  V  =  Vv  so  that 
^  . 

V  =  V  {cos $  [v  —  (i-s)s]  —  sin <j>vx.s  +  («•«)!}  .  (A12) 

For  a  given  v,  any  choice  of  s  —  provided  it  is  not  collinear  with  v  —  determines  a  right-handed 
orthogonal  Pauli-space  basis  q,f,s  in  the  following  way.  Take  f  to  be  the  direction  of  sxv,  i.e., 
sxv  =  sin0f  with  6  =  l(s,v).  Then,  since  v  is  orthogonal  to  f,  it  can  be  written  as 

v  =  sin  0  q  +  cos  6  s  , 

where  q  =  rxs.  Thus  Eq.  (A12)  becomes 

V'  =  V  (A13a) 

—  sin  0  (cos  <j)  q  -{•  sin  cos  0  s  ,  (Al3b) 

which  means  that  V  is  just  a  rotated  version  of  V .  In  standard  spherical  coordinates  relative  to 
the  q,f,s  axes,  v'  is  a  unit  vector  with  spherical  angles  (6,<f>),  whereas  v  was  a  unit  vector  with 
spherical  angles  (0,0).  Clearly,  since  V0'  =  Vb  and  V'  =  V,  the  entire  V  — ►  V'  transformation 
amounts  to  nothing  more  than  a  Pauli-space  rotation  of  the  vector  part  of  V  about  the  direction  s 
through  an  angle  <f>  —  —2 S.  Since  S  is  Hermitian,  S\,Si,  S3  are  all  real,  making  the  rotation  angle 
real  also. 

If  S  is  not  Hermitian,  then  the  magnitude  S,  as  defined  above,  can  become  imaginary,  giving 
the  unit  vector  s  both  real  and  imaginary  components.  A  more  general  class  of  “pseudo-Hermitian” 
S  operators  (with  real  S3  but  imaginary  Si,  S2 )  is  encountered  below.  The  above  result  still  applies, 
provided  the  rotation  angle  is  allowed  to  be  imaginary. 
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Appendix  B 

RELATION  TO  THE  BORN  SERIES 


This  appendix  relates  the  “Foldy-Wouthuysen”  (FW)  solution  obtained  in  the  body  of  the 
report  to  the  solution  achieved  in  a  scattering  context  by  the  Born  method.  The  two  are  shown  to 
agree  through  second  order. 

Scattering  Problem  in  Integral  Form 

The  Helmholtz  equation  of  interest  is 

dxii(x)  +  u(x)  =  —2 e(x)u(x)  ,  (Bl) 

in  the  scaled  x  coordinate.  A  scattering  problem  is  being  considered,  so  there  are  no  sources  at 
finite  x.  The  retarded  Green’s  function  for  a  uniform  medium  (one  with  e  =  0)  is 

g{x,y)=l-ei\*-y\  .  (B2) 

It  satisfies 

<Pxg{x,y)  + g{x,y)  =  -6(x  -  y)  (B3) 

and  allows  Eq.  (Bl)  to  be  re-expressed,  within  any  interval  a  <  x  <  b,  as 

u(x)  =  [g(x,  y)dyu(y)  -  u(y)dyg(x,y)]yy=a  +  /  dy2g(x,y)e(y)u(y)  .  (B4) 


Bounded  Scattering  Region 


For  simplicity,  suppose  that  all  environmental  nonuniformity  lies  in  a  bounded  “scattering 
region”  —  specifically  that  the  support  of  e(x)  is  the  finite  interval  l  <  x  <  r.  Outside  that  interval 
the  field  must  have  the  form 


(  Ae+i^x~a^  +  Be ...  x  <  i 
\  Ce+’t1"0*  +  De-l'(x-a)  ...  r  <  x 


(B5) 


where  the  point  x  =  a  has  been  used  for  the  (arbitrary)  phase  reference  and  A,  B ,  C,  D  are  param¬ 
eters.  One  can  always  arrange  for  A  =  1  and  D  =  0  by  normalizing  appropriately  and  taking  the 
incident  signal  to  come  from  the  left.  Then  if  a  <  l  and  r  <  b,  the  [•••]*  term  in  Eq.  (B4)  reduces 
to  e*(x~a),  leaving 

fb 

u(x)  =  et(x~a^+i  dye^x~v^e(y)u(y)  ...  a<x<b.  (B6) 
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Now  we  rename  a  =  x o  and  let  b  — *  oo  so  that  the  integral  equation  to  be  solved  is 

u( x)  =  uq(x )  +  i  f  dye^x~y^e(y)u(y)  ...  xo  <  x  <  oo  ,  (B7) 

Jx  0 

given  the  incident  signal  uq(x )  =  e^x~x°\ 


Born  Method 

The  smallness  of  e(x)  suggests  that  one  attempt  to  obtain  the  solution  u  as  a  Born  series 

=  u>o(a:)  +  u>i(z)  +  -  (B8) 

with  wn  =  0(e"),  starting  with  the  incident  field  wq  =  uq.  Born’s  method  obtains  the  partial  sums, 
un  =  wq  +  •  •  •  +  wn,  of  this  series  by  iteration, 

un(x)  =  u0(x)  +  i  f  d.ye'\x~y\t(y)un-i{y)  ...  n>  0.  (B9) 

In  other  words,  successive  terms  in  the  Born  series,  Eq.  (B8),  are  generated  via 

wn(x)  =  i  f  dy  e^x~v^e(y)wn-i(y)  ...  n  >  0  ,  (BIO) 

Jx  0 

and  the  nth  Born  approximation  is  uB  =  un  +  0(en+1). 

Preliminaries 


Before  proceeding  in  that  direction,  it  will  be  convenient  to  introduce  some  shorthand  notation 
and  a  preliminary  result. 


Differentiation  will  be  denoted  by  superscripts  in  parentheses, 


7(n)(x)  = 


dnq(x ) 
dxn 


...  n  >  0  , 


(Bll) 


which  is  a  departure  from  the  usage  in  the  body  of  the  report,  (n  =  0  is  the  trivial  case  =  q.) 
This  will  be  extended  to  n  <  0  so  that,  for  example,  q^1'1  denotes  the  anti-derivative, 


9(  1}(z)  =  /  dyq{y)  . 

Jxq 


In  this  notation,  a  typical  Taylor  series  takes  the  form 

00  „(n 
n! 


00  n 

q(x  +  y)  =  ZT<l{n)(x)  • 


(B12) 


(B13) 


n= 0 


We  will  also  need  the  following  result: 


n+1 


(B14) 


This  may  be  obtained  by  first  changing  the  integration  variable  to  t  =  2y,  thereby  converting  the 
left-hand  side  to 


1  /'OO 

-nrrr  /  • 

n!  2n+1  Jo 
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For  physical  reasons7,  that  singular  integral  is  to  be  understood  as  the  |-y |  — »  0  limit  of 

!°°  dte~ , 

JO 

which  can  be  found  in  tables  [Bl]: 

r&ir-'e- ll|,e“l  =  77^4«e’>arctan(</bl)  •••  »(/*)  >  -1,  \l\  >  |3(«)|  • 

Jo  (o^  +  Yr' 

With  y,  =  n  +  1  and  5  =  1,  the  I'yj  — ♦  0  limit  is  T(n  +  l)e’(n+1),r/2  =  n!  in+1,  confirming  Eq.  (B14). 
First  Born  Approximation 
The  first  Born  term  is 


roo 

Wl(x)  =  i  /  dyei{y-x°^x-y^e(y)  , 

Jx  o 


or,  with  the  integration  range  partitioned  at  y  =  x, 


wi(x)  = 


i  [  dye(y) 
Jxq 


e+i(x-x0)  _j_ 


w; 


dye 


»2(j/-x0) 


e(y) 


,-i(x-x  o) 


(B15) 


/(*) 


6(x) 


Outside  the  scattering  region,  where  the  functions  /  and  b  are  constant,  the  significance  of  the 
two  terms  in  Eq.  (B15)  is  unambiguous.  To  the  right  (for  r  <  x),  b  vanishes  and  the  first  term 
reduces  to  f(+oo)e+,^x~x°^  where 

f°C 

f(+oo)  =  i  dye(y)  , 

J  —  OC 

while  to  the  left  (for  x  <  £),  f  vanishes  and  the  second  term  reduces  to  &(— o o)e~^x~x°^  where 

/OO 

dy  el^y~x°^e(y)  . 

-OO 

On  the  right  side  of  the  scattering  region,  «q  is  a  constant-amplitude  wave  that  moves  to  the  right; 
on  the  left,  it  has  a  different  constant  amplitude  and  moves  to  the  left.  Note  also  that  for  x  <  £, 

/OO 

dy  ei2ye(y)  . 

'OO 

In  terms  of  the  unsealed  lengths  X  =  x/ko,  Y  =  y/ko,  this  is 

uiPO  =  e-ikoXo{eikoX  +  [ik0e(2k0)]e-ikoX}  , 

in  which  the  [•  •  •]  factor  is  recognizable  as  the  first-order  Born  reflection  coefficient  expressed  in  the 
standard  way  [B2]  as  a  Fourier  transform  e(/c)  =  dY  e'nY e(Y). 

7j7|  embodies  the  effect  of  attenuation  in  the  medium. 
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But  we  are  interested  mainly  in  the  field  within  the  scattering  region,  where  w\  does  not  yield  to 
any  such  simple,  unambiguous  reduction  into  forward  and  backward  components.  In  the  scattering 
region,  however,  Eq.  (B15)  can  be  written 


rx  roo 

wi(x)  =  wo(x)  i  /  dye(y)  +  i  /  dyet2ye(x  +  y)  . 
.  Jx o  Jo 


(B16) 


Then,  with  e(x  +  y)  =  [e(x  +  y)  —  e(x)]  +  e(x)  in  the  second  integrand  and  with  f//°  dye'2y  =  i/2 
(Eq.  (B14)  with  n  —  0),  one  has 


w\(x)  =  wq(x)  I'll  dy  e(y)  -  ^e(x) +  i  JQ  dV  e,2y[c(a:  +  y)  -  e(x)] 


(B17) 


This  form  emphasizes  the  fact  that  w\  can  always  be  regarded,  without  any  approximation  at  all, 
as  a  forward-propagating  carrier  wq(x)  =  exp{i(x  —  xo)}  with  an  x-dependent  modulation  y  +  u. 
It  also  suggests  that,  although  y  ~  e,  it  may  be  possible,  when  e  is  a  slowly  varying  function,  to 
relegate  v  to  the  0(e2)  terms  so  that  it  contributes  nothing  to  the  Born  approximation  at  this 
order.  The  remainder  of  this  appendix  is  basically  a  systematic  pursuit  of  this  possibility  through 
first  and  second  order. 


When  e(x  +  y)  in  the  second  integrand  in  Eq.  (B16)  is  expanded  in  a  Taylor  series  about  y  —  0, 
and  Eq.  (B14)  is  used  to  evaluate  the  resulting  /0°°  dy  el2yyn/n\  terms,  the  outcome  is 


oo  /  ,•  \  n+l 


wi/wq  =  i  jr  Q) 


(B18) 


(Note  that  y  comes  from  n  =  —1, 0  and  v  from  n  >  1.)  Thus  the  first  Born  partial  sum  is  given  by 


OO  /  ,•  \  n+l 
Ul/W0  =  1  +  i  Yl,  [2 )  6  " 


(B19) 


rx  1  00  / »'  \  fl+1 

(x)  =  e,(l_Io)  1  +  ij  dye(y)- 2e(a;)  +  iZj  (2)  e(n)(*) 


With 


^  ~  e 


the  first  Bom  approximation  uB  =  u\  +  0(e2)  can  then  be  written 


(B20) 


(B21) 


uB(x)  =  e*x-x^  ^1  -  ^e(x)^  e^o dyt(y)  +ijt  (i)”+1  €(")(x^  +  >  (B22) 

but  that  is  as  far  as  one  can  go  without  some  knowledge  of  the  behavior  of  the  derivative  terms  in 
the  remaining  sum.  If,  as  in  section  7,  one  assumes  that 


e(n)~en+1  ...  n  >  0  , 


(B23) 
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then  all  of  them  belong  with  the  0(e2)  terms  and  one  has 

«»(*)=  (l  -  !<(*))  e  £ * +  0(e2)  .  (B24) 

Term-by-term  comparison  with  the  first-order  form  of  Eq.  (110),  namely, 

uFW(x )  =  Q[1](x)e-«^1I(*)e<V[1](x,xo)  +  0(e 2)  ^  (B25) 

confirms  that  they  agree  through  first  order, 

uB  -  uFW  =  0(e2)  .  (B26) 

Second  Born  Approximation 


The  second  Born  series  term  is  given  by  the  double  integral 


roo  r  oo 

W2(x)/w o(x)  =  -  /  dy  dzE(y,z ) 

Jx  0  Jxq 

(B27) 

with  the  integrand8 

E(y,z)  =  e^y-x\+\z-y^z-x^e{y)e{z)  . 

(B28) 

When  the  integration  ranges  are  partitioned  at  y  =  x  and  z  =  y,  this  becomes 

w2/wq  =  —(a  +  b  +  c  +  d) 

(B29) 

\ 

in  terms  of  the  four  quantities 

rx  r  oo 

a(x)  =  dy  dz  e,2^z~y^e(y)e(z) 

Jx  o  Jy 

(B30a) 

b(x)=  [  dy  f  dz  e(y)c(z) 

Jx 0  Jx 0 

(B30b) 

/•oo  ry 

c(x)  =  /  dy  dz  e,2(y~x^e(y)e(z) 

Jx  Jx  0 

(B30c) 

roo  roo 

d(x)=  dy  dz  el2<'z~xh{y)e(z) 

Jx  Jy 

(B30d) 

which  will  now  be  evaluated  in  turn. 

The  a  term  is 

rx  roo 

a(x)  =  /  dyt{y)  /  dzel2ze(y  +  z) 

Jx  o  JO 

(B31a) 

rx  °°  roo  7n' 

=  /  dye(y)Y,t{n\y)  /  dze'2z  — 

Jx0  £^0  L/o  n\  \ 

(B31b) 

00  /,\n+l 

=  £  u)  <*"(*) 

n=0 

(B31c) 

8£’s  x  dependence  is  left  implicit. 
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where  expanding  e(y  4-  z)  in  a  Taylor  series  about  z  —  0  produces  Eq.  (B31b),  and  Eq.  (B31c) 
comes  from  evaluating  [•  •  •]  via  Eq.  (B14)  and  using  the  definition 


an(z)  =  [  dye(y)e^n\y)  . 
Jx  0 


(B32) 


The  b  term  is,  with  no  approximation, 

b (*)  -\\jx  dy€(y ) 


(B33) 


The  c  term  is 


where  the  definition 


:(*)  =  [  dye,2(-y  x)^{y) 

Jx 

=  J  dyet2y^{x  +  y) 

Lio  n! 


°°  /  •  \  n+1 

=  E(i)  W*) 

n=0 


j(x)  =  e(x)  [  dye(y) 

Jxq 


(B34a) 

(B34b) 

(B34c) 

(B34d) 

(B35) 


is  used  in  Eq.  (B34a),  then  expanding  j(x  +  y)  about  y  =  0  produces  Eq.  (B34c),  and  evaluation 
of  [•  •  •]  through  Eq.  (B14)  produces  Eq.  (B34d). 


The  d  term  is 


roo 

d(x)  —  /  dy  et2^y~x^a^\y) 

Jx 

f°° 

=  /  dy  e,2y  a^\x  +  y) 

JO 

°°  °°  /i\n+i  ,  .  r  /■=»  7/i 

-ssG)  L 

=  EEG)n+£+24m)(-) 

n=0£=0  V  ' 
oo  /  «  \  n+2  n 

-E(j)  Z°l+e%) 

n=0  v  '  e=o 


(B36a) 

(B36b) 

(B36c) 

(B36d) 

(B36e) 


where  expanding  a^\x  +  y)  about  y  =  0  and  using  the  definition  Eq.  (B32)  produce  Eq.  (B36c), 
then  evaluating  [•  •  •]  using  Eq.  (B14)  produces  Eq.  (B36d),  and  summation  in  an  alternate  order 
yields  Eq.  (B36e). 
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The  second  Born  partial  sum  is  given  by 


112/wo  =  1  +  wi/wo  +  W2/W0 


(B37) 


with  wi/vjq  and  W2/wq  obtained  from  Eq.  (B18)  and  Eq.  (B29),  respectively.  For  the  second  Born 
approximation,  we  will  need  to  identify  and  retain  contributions  through  order  e2.  For  wi/wq  that 
is  easy;  to  second  order  it  is  simply 

4 

For  W2/wq,  we  need  all  the  second-order  terms  from  a,  b,  c,  and  d.  Since  b  is  purely  second-order, 
its  contribution  is  simply  e(-1)2/ 2. 

For  a,  we  need  to  assess  the  an  =  ^  terms,  the  first  two  of  which  are 

oto(x)=  [  dye2(y ) 

Jx  0 

ai(x)=  [  dye(y)e(x\y)=]-e2{y) 

Jxq  &  Xh 


Since  e(")(xo)  =  0  for  all  n  >  0,  we  have 


ao  =  ^  ^  ~  e2 
<*1  =  ^e2  ~  e2  , 


(B38) 

(B39) 


both  of  which  contribute  at  order  e2.  The  an  for  n  =  2, 3,  •  •  •  can  be  evaluated  by  repeated 
integration  by  parts.  For  even  and  odd  n,  respectively,  they  are 

a2e(x)  =  g(-)V*)(*)«("-*"1)(«)  +  (-)<  fXdye^2(y) 

tl  _ ,  _ _ ^ - ' 

^€2<+l 

««+!(*)  =  E(-)^(fc)(x)e(2^(x)  +  (-)^>2(x) 


for  i  >  1.  Since  2t  +  1  >  3,  all  of  these  belong  in  0(e3)  and  we  are  left  with 

a  =  ^<*0  -  +  0(e3) 


=  r2(_1)-r2+°(e") 


(B40) 


For  c,  we  need  the  7^  =  B  ^  ^  terms.  Since 


7(°)  _  e(0)e(-i) 


(B41) 
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these  can  easily  be  obtained  from  the  standard  expression  for  the  nth  derivative  of  a  product, 

7(n)  =  E(  )eWe(n_fc_1)  •••  n>0- 


With  the  k  =  n  term  written  separately,  this  is 


7(")  =  e(fc)*(n-*_1)  • .  •  n  >  0  . 

^,n+2  k=0  \  / 


(B42) 


fn+l 


Clearly,  the  only  contribution  below  the  third  order  is  the  k  =  0_one  for  n  =  1.  Thus 

c  =  ^70  -  ^71  +  0(e3) 

=  le(0)€(-l)_  l£(0)2+O(£3}' 


(B43) 


written  as 


The  fact  that  ^  =  |c(°)c(m)J^  ^  is  the  fth  derivative  of  a  product  allows  it  to  be 


|  e(k)e(m+e-k) 


Contributions  to  d  come  from  m  =  n  —  £.  Since  these  are 

J*+1)  _  Vs  l  1  \  „(*)>-*) 


<-t  =  E  it  ...  n>0, 


k= 0 


«»+2 


only  n  =  0  contributes  below  the  third  order.  Thus 

.  2 


«*=(“)  41>  +  0(e3) 

=  _Ie(0)2_J_o(e3) 

4 


(B44) 


(B45) 


In  light  of  the  above  results,  the  second  Born  approximation  is  given  by 

«BM,  =  1  +  -  ^2<-1)  +  §<*  -  +0(£3)  , 


(B46) 


or,  in  the  conventional  notation, 
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Term-by-term  comparison  with  the  second-order  form  of  Eq.  (110),  i.e., 

uFW(x)  =  QW(x)e-idl2]Weivl2](~x'Xo)  +  0(e3) 

=  (l  -  i£ (*)  +  |c2(x))  6 +  0(e3)  (B48) 

confirms  their  agreement  through  second  order, 

uB  -  uFW  =  0(e3)  .  (B49) 

To  extend  this  investigation  to  the  third  Born  approximation  would  mean  analyzing  w^/wq  — 
a  triple  integral  analogous  to  Eq.  (B27).  The  analog  of  Eq.  (B30)  would  involve  23  =  8  terms; 
furthermore,  their  third-order  contributions,  along  with  those  of  w^/wq  and  w\/wq,  would  need  to 
be  identified  and  retained.  The  effort  appears  prohibitive,  but  it  seems  compelling  to  conjecture 
that  uB  —  uFW  vanishes  to  all  orders. 
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